Phonon Hall thermal conductivity from Green-Kubo formula 
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We derive a formula for the thermal conductivity tensor of a ballistic phonon Hall model. It is 
found that, although the diagonal elements of the conductivity tensor diverge to infinite, the off- 
diagonal elements are finite, antisymmetric, and odd in magnetic field. The off-diagonal elements 
are non-zero only if the dynamic matrix of the phonon system breaks mirror refiection symmetry. 
The results are obtained without perturbative assumptions about the spin-phonon interactions. 
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The Hall effect of electronic conduction is well-known 
and has many applications. The analogous effects for the 
transport of gas molecules, spins, and photons also exit 
[Hi S [3- The phonon Hall effect, that is, the appear- 
ance of a transverse thermal current when a magnetic 
field is applied perpendicular to the direction of tem- 
perature gradient, is esoteric and not well understood. 
Electrons couple directly to the magnetic field through 
the Lorentz force. There is no obvious coupling between 
phonons and magnetic field. In 2005, Strohm et al. re- 
ported such an effect in a paramagnetic dielectric garnet 
Tb3Ga50i2 t^] and confirmed also in Ref. [5], and called 
it phonon Hall effect. Two theoretical papers followed 
@) Both of them considered a similar model of the 
spin-phonon interaction, and both of them treated the 
interaction perturbatively. The work of Ref. [3] appears 
to imply that ballistic systems cannot produce a phonon 
Hall effect, and the authors evoked further high order 
spin-phonon interaction terms to demonstrate the exis- 
tence of the effect. Although the two approaches are 
quite different, one based on Green-Kubo formula, the 
other on Boltzmann-type kinetic equation, curiously, the 
final results for the off-diagonal thermal conductivity ten- 
sor are similar. 

In this paper, we address the following issues. (1) Is 
a ballistic system capable of producing the phonon Hall 
effect? Our answer to this question is affirmative. Al- 
though the effect will be smaller as the linear size L of 
the system becomes larger (scaled as l/L). (2) What is 
the role of symmetry? We found that break of a mirror 
reflection symmetry is essential to observe the phonon 
Hall effect. If the system looks the same inside a mir- 
ror, we should not observe such effect on very general 
ground. We use the same model as that of Refs. @, Q 
but without the perturbative assumption. The pertur- 
bative expansion with respect to the spin-phonon inter- 
action breaks down near the F-point of the phonon dis- 
persion. This complicates the behavior of thermal con- 
ductivity at very low temperatures. Since the model is 
ballistic, the thermal conductivity in general should di- 
verge with the system sizes. But for isotropic systems like 
the two-dimensional square or honeycomb lattices, the 



off-diagonal thermal conductivity is in fact finite. In the 
following, we introduce the model, outline a derivation 
of the thermal conductivity using Green-Kubo formula, 
and present numerical results and give some comparison 
with experiments. 

We consider a harmonic periodic lattice with the ex- 
tra Rahman (or spin-orbit) interaction at each lattice site 
proportional to s- (r x p). Here s is the (pseudo-) spin rep- 
resenting the Kramer doublet; r and p are displacement 
and conjugate momentum. We'll replace s by an average 
magnetization of the system and choose the vector to be 
in z direction. The explicit spin degrees of freedom drop 
out of the problem. The Hamiltonian of the system can 
be written in a compact form 



H 



1 



1 



:P P 
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where m is a column vector of displacements away from 
lattice equilibrium positions for all the degrees of free- 
doms, multiplied by \/m, p is the associated conjugate 
momenta. The Rahman term, u^Ap, is onsite; the ma- 
trix A is an antisymmetric real matrix, A'^ = —A, and 
is block diagonal with diagonal elements (in two dimen- 
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We'll call h magnetic field although h is only propor- 
tional to the real magnetic field in a paramagnet. It has 
the dimension of frequency. Since the interaction term 
depends on momentum, the velocity and momentum are 
not the same but related through ii = p — Au. This is the 
same model studied in Refs. [a, 0] except a slightly dif- 
ferent notation. It has been proposed (in a more general 
form) based on quantum theory and fundamental sym- 
metries long time ago to study spin-phonon interactions 

Equation ([1]) is quadratic in the dynamic variables u 
and p, thus is amenable for an exact solution. Our calcu- 
lation procedure is as follows. We first obtain the eigen 
modes of the system. Using the eigen modes, we give 
an expression for the energy current. We then apply the 
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Green-Kubo formula to compute the thermal conductiv- 
ity tensor. Since the system is periodic, we can apply the 
Bloch theorem. The polarization vector e then satisfies 



(3) 



where -D(k) = dynamic ma- 

trix. Kiji is the submatrix between unit cell I and /' in 
the full spring constant matrix K] R; is the real-space 
lattice vector. This equation is not a standard eigen- 
value problem. It is numerically more advantageous to 
consider both the coordinates and momenta and to solve 
an eigenvalue problem: 



A D 

-I A 



(4) 



where x = {ji, e)^ and / is an identity matrix. Contrary 
to the usual lattice dynamic problems, the polarization 
vectors are not orthogonal to each other. We need to 
consider both the right and left eigen vectors. Because of 
the special form of Eq. ([4]), the left eigen vectors and right 
eigen vectors are not really independent. It is possible to 
choose the left eigenvectors x = {fl,e) = [e^—^^). The 
orthonormal condition then holds between the left and 
right eigen vectors. In particular, the eigen modes can 
be normalized according to 
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Since the matrix on the right-hand side of Eq. ^ is not 
anti-hermitian, there is no guarantee that the frequencies 
uj will be real, but the eigenvalues always come in pairs, 
ztiv. We take only ui > modes. With these choices of 
the eigen modes, displacement and momentum operators 
can be taken in the standard second quantization form. 
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ak + h.c. 
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where k — (k, a) specifies the wavevector as well as 
phonon branch, ak is the annihilation operator, and h.c. 
stands for hermitian conjugate. The momentum and 
displacement polarization vectors are related by, e.g., 
fi = —iuoe + Ae. We can verify that the canonical com- 
mutation relations are satisfied, [ui,pj,] — ihSijiI, and 

H = Ek^Maick + 1/2). 

Based on a definition of the local energy density and 
the continuity equation for energy conservation, an en- 
ergy current density can be defined as @, 0, [HI , 
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where the index c = x, y, or z labels the cartesian axis, 
V is the total volume of N unit cells. The components of 
the current density vector can be expressed in terms of 
the creation/annihilation operators. Ignoring the a a and 
a^a^ terms which vary rapidly with time, one obtains Q 
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The thermal conductivity tensor can be obtained from 
the Green-Kubo formula , 
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where (3 = l/{kBT), the average is over the equilibrium 
ensemble with Hamiltonian H. The time dependence of 
the annihilation operator is trivially given by ak(t) = 
flfce""^''*. This is also true if t is imaginary. Substituting 
the expression J'^ into Eq. pO)) . using the result 
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where fi — (e'^'*'^' — 1) ^ is the Bose distribution function, 
we obtain 
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where the F function is defined as 
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To simplify notations, we have suppressed indices, e.g., 
UJ = a.'cr(k), e' — e(j'(k). We have added a damping 
term e"*** when integrating the oscillatory factor. The 
diagonal element of F is related to the group velocity, 
F^^(k) = 2dujl/dk''. The off-diagonal elements are m 
general not zero. The first term in Eq. (|lip factors into 
two independent summations which does not contribute 
to Kab due to symmetry of with respect to the 

wavevector k. Equation ((T^), together with the defini- 
tion (|13[) . is the main result of this paper. 

We make some general comments on Eq. p2|). The 
first and last factors inside the summation sign can be 
combined, (e'"*('^'-'^) - l)/'(/ + 1) = / - /'. Written in 
this way, the equation resembles the Landauer formula 
for ballistic transport. The second factor makes the con- 
ductivity diverge in the form 1 /r] unless the leading term 
in an expansion in the damping factor rj happens to be 
zero. The diagonal elements Kaa indeed diverge to infi- 
nite. This is expected, as the system is ballistic consisting 
of independent oscillating modes. There is no intrinsic 
scattering mechanism in the system. 
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FIG. 1: Phonon dispersion relation of a triangular lattice, (a) 
The angular frequency of longitudinal mode as a function of 
kya with fcj; = 0. The bottom curve is /i = 0; top curve is 
h — 5x 10^^ rad s~^ . (b) The frequency as a function of at a 
fixed wavevector ka — (0, 1). The top curve is the longitudinal 
mode, the bottom (broken) curve is the transverse mode. 

The ofF-diagonal elements do not diverge if the system 
is isotropic in the sense that Kab is independent of the 
choice of the coordinate axis. In this special case, the 
ofF-diagonal elements are antisymmetric and odd in the 
magnetic field h, Kabih) = —Ki,a{h) ~ K.ba{—h), consis- 
tent with the Onsager relation. This property does not 
hold for arbitrary anisotropic systems. We argue that in 
the isotropic case, Eq. (fT2|) is physical and is the correct 
prediction for the Hall thermal conductivity. 

Even in the isotropic case, the off-diagonal term is zero 
unless reflection symmetry is broken. More precisely, if 
there exists an orthogonal transformation independent of 
k such that SDS^ = D, SAS^ = -A, then Kab = for 
a b. The physical meaning of this symmetry is clear. If 
we look the system in a mirror, since D is the same and A 
flips a sign, but the physics must be invariant, we should 
have Kab{D,A) — Kab{D, —A). But Kab{D,A) must be 
an odd function in A. So we must have Kab = 0, a b. 
This property should be quite general, independent of the 
models used. As an example of systems with vanishing 
ofF-diagonal thermal conductivity, we can cite a square 
lattice (or cubic lattice) with only the nearest neighbor 
coupling with a dynamic matrix which is diagonal. 

In the following, we present numerical results based on 
Eq. p2|) . But first, we discuss some interesting features 
of the phonon dispersion when the Raman interaction 
term is turned on. In Fig. [l] we show results for a tri- 
angular lattice with only the nearest neighbor couplings. 
The coupling matrix between two sites is such that the 
longitudinal spring constant is — 0.144 eV/(uA^) and 
the transverse spring constant Kt is 4 times smaller. The 
unit cell lattice vectors are (a,0) and (a/2,a\/3/2) with 
a = 1 A. This choice gives longitudinal and transverse 
sound velocity 3981 m/s and 1921 m/s, respectively, com- 




h (10'"rad/s) 

FIG. 2: Thermal Hall conductivity as a function of the cou- 
pling h for fixed temperatures T — 5K and 100 K, respec- 
tively. 

parable to typical experimental values. At the F-point, 
the effect of the interaction is to shift the frequencies 
from Wo to ujo zt h, for both the acoustic modes and op- 
tical modes (if any). In particular, the acoustic modes 
develop a gap from 0. Away from the F-point, the cor- 
rections are of order h"^. Due to the interaction, some 
modes have imaginary frequencies and are no longer sta- 
ble. This is very pronounced for the transverse modes 
for large h, see Fig. 1(b). The system can be stablized, 
at least for small h, by adding a small onsite potential 
(which, of course, breaks the translational invariance of 
the lattice). If we change the model to use velocity ii 
instead of the conjugate momentum p in the interaction 
term, an onsite term of magnitude h'^ is generated natu- 
rally. In such a model, these instabilities do not appear. 
However, there is no good reason to use vFAii instead of 
u^Ap other than the above observation. 

In Fig. 2, we give the off-diagonal thermal conductivity 
Kxy of the triangular lattice (assuming 1 A thick) as a 
function of h for two different temperatures, T = 5 K 
and 100 K. It is clear that, for small h, the dependence 
of on h is linear. For large h, the growth becomes 
weaker than linear. For very large h (not shown), due 
to the instability, K^y becomes rather singular, and can 
even become negative. This range of parameters is not 
physical. In computing the results of Fig. 2, we have 
added a small onsite value of order 10~^i^j;,. The results 
are sensitive for this onsite value only for large h, but are 
nearly independent of the onsite value for small h. 

In Fig. 3, we display the temperature dependence of 
the off-diagonal thermal conductivity K^y It is seen that 
Kxy saturates at about 100 K a,t h = O.lradTHz. At low 
temperatures, Kxy decreases with temperature approxi- 
mately linearly. However, due to a complicated effect of 
h to the dispersion relation, it appears that Kxy has a 
faster fall off than linear. 

We comment on experimental data 0, Q in compar- 
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h= l.UradTHz 
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FIG. 3: Thermal Hall conductivity function of the 

temperature T for fixed coupling h = lO^'^rad/s (solid line) 
and lO^^rad/s (dotted line). 

ison with our numerical results. We have not taken 
into account the specific lattice structure and atomic 
details used in experiments. A quantitative compari- 
son is not possible. However, the phonon model pa- 
rameters are comparable to real systems by matching 
the sound velocities. The most uncertainty in a com- 
parison is the coupling h. The experimental value for 
K^y at T = 5.13 K and magnetic field H = 3T is 
2.0 X 10~^Wm~^K~^ 0. This is consistent with a very 
small value of ft, = 10~'^radTHz. Although the diago- 
nal element k^^ {— i^yy) diverges to infinite in our the- 
ory, we can choose a finite rj in Eq. (|12p to mimic a fi- 
nite phonon life time. We find a very weak dependence 
of Kxx on h. On the scale of ft ~ IradTHz, k^x is 
nearly a constant. By matching the experimental value 
of order 0.5 Wm~^K~^, we can infer a mean free path 
I — cjr] ~ 10^ A (where c is sound velocity), which ap- 
pears a bit too small given the very low temperatures in 
experiments. 

It is interesting to compare the present treatment with 
that of nonequilibrium Green's function (NEGF) ap- 
proach in Ref. [I3i| . The qualitative features are in agree- 
ment, such as the vanishing phonon Hall effect on square 
lattice. In NEGF approach, the leads are modeled ex- 
plicitly. It was assumed that leads do not have the spin- 
phonon interaction. This has the advantage of stablizing 
the system, even though the spin-phonon system repre- 
sented by the Hamiltonian if, Eq. ([1]), may be unstable 
as a bulk system. NEGF deals with very small systems 
in practice. Some of the oscillatory behavior, perhaps of 
a reflection of the wave nature, is not found here. The 
present theory is more suitable for comparison with ex- 
periments which were done on samples of mm scale. 

Another point is the role of nonlinear interactions. The 
phonon-phonon and spin-phonon interactions will pro- 
duce a finite life time for the phonons, rendering a finite 



thermal conductivity tensor for all components. We ex- 
pect that, if there is a systematic expansion in terms 
of the phonon life-time or in terms of the interaction 
strength, our main result, Eq. (fT2|) . should be the lead- 
ing contribution. The interaction should give only small 
corrections. 

In summary, we have presented a theory of phonon 
Hall effect based on a ballistic lattice dynamic model. 
It is shown that the phonon Hall effect can be present, 
provided that the system does not possess a reflection 
symmetry. This is different from Ref. , which does not 
suggest this subtlety. Since the Hamiltonian is quadratic 
in the dynamic variables, a perturbative treatment is not 
necessary. In fact, it fails near the F-point. We have 
given numerical results on a simple two-dimensional tri- 
angular lattice and the qualitative features are the same 
for all lattices in two and three dimensions. In partic- 
ular, it is not necessary that the system has degenerate 
phonon modes when ft = 0. When more elaborate model 
is known (e.g., from a first-principles calculation), the 
current theory can be applied to more realistic systems. 
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